% load parameters (posterior mean) 


alpha_1    =      1.79599541;  
alpha_2    =      1.14322405;  
gamma_1    =      0.09410810;  
gamma_2    =      0.08092049;  
beta       =      0.99914173;  
kappa      =      0.44804350;  
tau        =      3.98004910;  
rho_a      =      0.04901146;  
rho_f      =      0.94985429;  
rho_i      =      0.68761632;  
sigma_a1   =      0.01528579;  
sigma_f1   =      0.00325211;  
sigma_i1   =      0.00614372;  
sigma_a2   =      0.00693212;  
sigma_f2   =      0.00115753;  
sigma_i2   =      0.00324415;  
pi         =      0.00723860;  
gy         =      0.00496256; 
p_11       =      0.96664397;  
p_22       =      0.94965365;  
q_11       =      0.90547855;  
q_22       =      0.93909420;  

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------
alpha1_grid = 0.85:0.05:3; 

mhrund    = 1; 


%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec2;

% parameter names
pnames = feval(fnames_.pnames,msel);


%--------------------------------------------------------------------
% means and standard deviations
%--------------------------------------------------------------------

% load file
	lfile = [path_.para 'stat_mhrun1'];
	load(lfile);

%p22_grid = 0.5:0.01:0.99; 
%gammagrid=0.01:0.01:0.5;

lk = length(alpha1_grid); 
para = drawmean';
wu1 = zeros(lk,1); 
wu2 = zeros(lk,1); 
wu3 = zeros(lk,1); 
wu4 = zeros(lk,1);
infh =zeros(lk,1);
infl = zeros(lk,1); 
for j=1:lk 
        
alpha1  = alpha1_grid(j);
%alpha1 = 1.79599541; 
alpha2 = 1.14322405; 
%gamma1 = gammagrid(j);
gamma1 = 0.09410810; 
gamma2 = 0.08092049; 
beta1  = 0.99914173; 
kappa1 = 0.44804350; 
tau1   = 3.98004910; 
rho_a  = 0.04901146; 
rho_f  = 0.94985429; 
rho_i  = 0.68761632; 
sd_a1  = 0.01528579; 
sd_f1  = 0.00325211; 
sd_i1  = 0.00614372; 
sd_a2  = 0.00693212; 
sd_f2  = 0.00115753; 
sd_i2  = 0.00324415; 
piess  = 0.00723860; 
gy     = 0.00496256; 
lnA_0  = 9.52881130; 
yst    = 0.06727340; 
%p11    = p11_grid(j); 
p11    = 0.96664397; 
%p22    = p22_grid(j);
p22    = 0.94965365; 
p11_f  = 0.90547855; 
p22_f  = 0.93909420; 




P_i = [p11 1-p11;
       1-p22 p22];
   
P_f = [p11_f 1-p11_f;
       1-p22_f p22_f];   

P = kron(P_i,P_f);


%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------

  A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
        -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];

    solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

    A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
        -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];
    solvec_f = inv(A_Af_mat) * [0;0;-kappa1;-kappa1;0;0];

    A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
        -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];
    solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];

   
    pi1a = solvec_a(3);
    pi1f = solvec_f(3);
    pi1i = solvec_i(3);
   

    var_a1 = (sd_a1^2)/(1-rho_a^2);
    var_u1 = (sd_f1^2)/(1-rho_f^2);
    var_e1 = (sd_i1^2)/(1-rho_i^2);

    var_a2 = (sd_a2^2)/(1-rho_a^2);
    var_u2 = (sd_f2^2)/(1-rho_f^2);
    var_e2 = (sd_i2^2)/(1-rho_i^2);
    
    wu1(j)  =((pi1f^2)*var_u1)/(((pi1a^2)*var_a1)+((pi1f^2)*var_u1)+((pi1i^2)*var_e1));

    wu2(j)  =((pi1i^2)*var_e1)/(((pi1a^2)*var_a1)+((pi1f^2)*var_u1)+((pi1i^2)*var_e1));
   
    wu3(j)  =((pi1f^2)*var_u2)/(((pi1a^2)*var_a2)+((pi1f^2)*var_u2)+((pi1i^2)*var_e2));  
    
    wu4(j)  =((pi1i^2)*var_e2)/(((pi1a^2)*var_a2)+((pi1f^2)*var_u2)+((pi1i^2)*var_e2));  
   
    infh(j) = wu1(j)*rho_f+wu2(j)*rho_i+(1-wu1(j)-wu2(j))*rho_a;
    
    infl(j) = wu3(j)*rho_f+wu4(j)*rho_i+(1-wu3(j)-wu4(j))*rho_a;
    
end


linewidth=2;
linecolor1='r';
linestyle1='-';
linestyle2='--';
%subplot(2,1,1)
xxtp1=   1.79599541;  
xxtp2=   1.14322405; 
ymin1=0.7; ymax1=infh(20);
ymin2=0.7; ymax2=infh(7); 
plot(alpha1_grid,[infh infl],'k')
ylim([0.7 0.95]); 
line([xxtp2 xxtp2],[ymin2 ymax2],'LineStyle',linestyle2,'Color',linecolor1,'LineWidth',1);
line([xxtp1 xxtp1],[ymin1 ymax1],'LineStyle',linestyle1,'Color',linecolor1,'LineWidth',1);
% subplot(2,1,2)
% plot(gammagrid,wu2,'k')
% 
% line([xxtp2 xxtp2],[ymin2 ymax2],'LineStyle',linestyle2,'Color',linecolor1,'LineWidth',1);
% line([xxtp1 xxtp1],[ymin2 ymax2],'LineStyle',linestyle1,'Color',linecolor1,'LineWidth',1);

			